MODULE  SUCOTWO_MOD
CONTAINS
SUBROUTINE SUCOTWO(PRVR0VT,YDVEG,YDCST,YDAGS)
!***
! (C) Copyright 1997- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.


!**   *SUCOTWO* - DOES THE INITIALISATION OF AGS PARAMETERS

!     A. Boone       * Meteo-France *     27/10/97
!     (following Belair)
!     MODIFIED BY
!     V. Rivalland                        10/04/02	 
!        Add: PAH and PBH coefficients for herbaceous water stress response
!     P. LeMoigne                            03/04   
!        computation of zgmest in SI units
!     M.H. Voogt (KNMI) "C-Tessel"  09/2005 
!      S. Lafont (ECMWF) externalised CTESSEL 04/2006 
!     Boussetta/Balsamo 11/2010 clean and compact code
!     Boussetta/Balsamo 07/2011 MACC budget respiration optimisation
!     M. Kelbling and S. Thober (UFZ) 11/6/2020 use of parameter values defined in namelist 
!     Agusti-Panareda 09/2021 Use optimized respiration from namelist
!     I. Ayan-Miguez (BSC) Oct 2023: move derived spatially distributed parameters to surf/module/susdp_deriv_ctl_mod.F90 routine

!    PURPOSE
!    -------

!     Initialize model to calculate net assimilation of 
!     CO2 and leaf conductance.
!              
!     METHOD
!    ------
!     Calvet at al (1998) [from model of Jacobs(1994)]

!     EXTERNAL
!     --------
!     none

!     REFERENCE
!     ---------

!     Calvet et al. (1998)

!     -------------------------------------------------------------------------
!
USE PARKIND1        ,ONLY : JPIM ,JPRB, JPRD
USE YOMHOOK         ,ONLY : LHOOK,   DR_HOOK, JPHOOK
USE YOS_VEG         ,ONLY : TVEG
USE YOS_CST         ,ONLY : TCST
USE YOS_AGS         ,ONLY : TAGS
USE COTWO_MOD

IMPLICIT NONE

REAL(KIND=JPRB)   , INTENT(IN)    :: PRVR0VT(:,:)
TYPE(TVEG)        , INTENT(IN)    :: YDVEG
TYPE(TCST)        , INTENT(IN)    :: YDCST
TYPE(TAGS)        , INTENT(INOUT) :: YDAGS

!*         0.     LOCAL VARIABLES.
!                 ----- ----------


INTEGER(KIND=JPIM) :: JVT    ! indexes for loops
INTEGER(KIND=JPIM) :: ICO2TYP  ! type of CO2 vegetation
INTEGER(KIND=JPIM) :: ZIDIA,ZFDIA,ZLON
INTEGER(KIND=JPIM) :: JI, JJ, JK    ! loop indexes
INTEGER(KIND=JPIM) :: IM, KN   
LOGICAL            :: ZLDLAND(1)                               
REAL(KIND=JPRB)    :: ZCO2MWG,ZAIRMWG
REAL(KIND=JPRB)    :: ZX1
REAL(KIND=JPRB)    :: ZX2
REAL(KIND=JPRB)    :: ZXM, ZXL, ZZ, ZP1, ZP2, ZP3, ZPP, ZZ1   
REAL(KIND=JPRB), PARAMETER   :: ZPEPS = 3.0E-14_JPRB  ! convergence criteria
     
   !    ----------------------------------------------------------------

! as a function of photosynthesis mechanism, C3=>1, C4=>2
REAL(KIND=JPRB), PARAMETER, DIMENSION(2) :: &
 & ZTOPT    = (/25.0, 35.0/), &
 & ZFZERO   = (/0.85, 0.50/), & 
 & ZFZEROST = (/0.95, 0.60/), &
 & ZEPSO    = (/1.417E-8, 1.167E-8/), & 
 & ZGAMM    = (/6.851E-5, 4.263E-6/), & 
 & ZQDGAMM  = (/1.5, 1.5/), & 
 & ZQDGMES  = (/2.0, 2.0/), & 
 & ZT1GMES  = (/5.0, 13.0/), & 
 & ZT2GMES  = (/36.0, 36.0/), &  
 & ZAMMAX   = (/1.833E-6, 1.417E-6/), & 
 & ZQDAMMAX = (/2.0, 2.0/), & 
 & ZT1AMMAX = (/8.0, 13.0/), &  
 & ZT2AMMAX = (/38.0, 38.0/), & 
 & ZAH      = (/2.381, 5.323/), &  ! Calvet (2000)
 & ZBH      = (/-0.6103, -0.8929/) ! Calvet (2000)

REAL(KIND=JPHOOK) :: ZHOOK_HANDLE
!    --------------------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SUCOTWO_MOD:SUCOTWO',0,ZHOOK_HANDLE)
ASSOCIATE(RABC=>YDAGS%RABC, RAIRTOH2O=>YDAGS%RAIRTOH2O, RAMMIN=>YDAGS%RAMMIN, &
 & RCO2=>YDAGS%RCO2, RCO2TOH2O=>YDAGS%RCO2TOH2O, RDMAX_AGS=>YDAGS%RDMAX_AGS, &
 & RMAIR=>YDAGS%RMAIR, RMC=>YDAGS%RMC, RMCO2=>YDAGS%RMCO2, RPOI=>YDAGS%RPOI, &
 & RRACCF=>YDAGS%RRACCF, RXBOMEGA=>YDAGS%RXBOMEGA, RCO2FRAC=>YDAGS%RCO2FRAC, &
 & RDAY=>YDCST%RDAY, RPI=>YDCST%RPI, &
 & NVTYPES=>YDVEG%NVTYPES)

!*         1.0 INITIALIZE local variables

ZXM=0.0_JPRB
ZXL=0.0_JPRB
ZZ=0.0_JPRB
ZP1=0.0_JPRB
ZP2=0.0_JPRB
ZP3=0.0_JPRB
ZPP=0.0_JPRB
ZZ1=0.0_JPRB

ZIDIA=1_JPIM
ZFDIA=1_JPIM
ZLON=1_JPIM
ZLDLAND=.TRUE.

!*         1.1 INITIALIZE CONSTANTS 
!              ---------- ---------

RMAIR = 28.970E-3_JPRD
RMCO2 = 44.011E-3_JPRD
RMC = 12.0E-3_JPRD

! atmospheric CO2 concentration (in kgCO2 kgAir-1) from IFS .
!  
ZAIRMWG = 28.970_JPRB
ZCO2MWG = 44.011_JPRB

!orig RCO2 = 0.000353_JPRB*ZCO2MWG/ZAIRMWG
!RCO2 = 0.000400_JPRB*ZCO2MWG/ZAIRMWG

RCO2 = RCO2FRAC*ZCO2MWG/ZAIRMWG

RDMAX_AGS = 0.045_JPRB
RRACCF = 1.00_JPRB
! if leaf single scattering albedo 'omega'=0.2, the foliage scattering
!      coefficient is

RAMMIN = 1.E-10_JPRB
RAIRTOH2O = 1.606_JPRB
RCO2TOH2O = 1.6_JPRB

! CH4 emission scaling factor depending on latitudinal band (North/Tropics/South)
! optimized at LSCE
! reoptimization at LSCE done on the final stage of CAMS41 (July 2021)
YDAGS%RVCH4S(1) = 0.00422_JPRB
YDAGS%RVCH4S(2) = 0.00301_JPRB
YDAGS%RVCH4S(3) = 0.00326_JPRB

!    -------------------------------------------------------------------------

!*         1.2 INITIALIZE PARAMETERS DEPENDING ON VEGETATION TYPE
!              ---------- ---------- --------- -- ---------- ----

IF (.NOT.ALLOCATED(YDAGS%LVSTRESS)) ALLOCATE (YDAGS%LVSTRESS(0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVBSLAI)) ALLOCATE (YDAGS%RVBSLAI(0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVLAIMIN)) ALLOCATE (YDAGS%RVLAIMIN(0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVSEFOLD)) ALLOCATE (YDAGS%RVSEFOLD(0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVGMES)) ALLOCATE (YDAGS%RVGMES(0:NVTYPES))
IF (.NOT.ALLOCATED(YDAGS%RVGC)) ALLOCATE (YDAGS%RVGC(0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVDMAX)) ALLOCATE (YDAGS%RVDMAX  (0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVF2I)) ALLOCATE (YDAGS%RVF2I(0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVCE)) ALLOCATE (YDAGS%RVCE(0:NVTYPES))
IF (.NOT.ALLOCATED(YDAGS%RVCF)) ALLOCATE (YDAGS%RVCF(0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVCNA)) ALLOCATE (YDAGS%RVCNA(0:NVTYPES))
IF (.NOT.ALLOCATED(YDAGS%RVANMAX)) ALLOCATE (YDAGS%RVANMAX(0:NVTYPES))
IF (.NOT.ALLOCATED(YDAGS%RXBOMEGAM)) ALLOCATE (YDAGS%RXBOMEGAM(0:NVTYPES))
IF (.NOT.ALLOCATED(YDAGS%RVR0VT)) ALLOCATE (YDAGS%RVR0VT(0:NVTYPES,2))
IF (.NOT.ALLOCATED(YDAGS%RVCH4QVT)) ALLOCATE (YDAGS%RVCH4QVT(0:NVTYPES))

! BATS table matching with ECOCLIMAP table 
! (1)  ! Crops, Mixed Farming			=>! 6  C3 CROPS
! (2)  ! Short Grass				    =>! 4  C3 GRASS
! (3)  ! Evergreen Needleleaf Trees		=>! 2  CONIFEROUS
! (4)  ! Deciduous Needleleaf Trees		=>! 2  CONIFEROUS
! (5)  ! Deciduous Broadleaf Trees		=>! 1  DECIDUOUS
! (6)  ! Evergreen Broadleaf Trees		=>! 3  EVERGREEN
! (7)  ! Tall Grass				        =>! 5  C4 GRASS7
! (8)  ! Desert					        =>! 
! (9)  ! Tundra					        =>! 4  C3 GRASS
! (10) ! Irrigated Crops			    =>! 6  C3 CROPS
! (11) ! Semidesert				        =>! 4  C3 GRASS
! (12) ! Ice Caps and Glaciers			=>!
! (13) ! Bogs and Marshes			    =>! 4  C3 GRASS
! (14) ! Inland Water				    =>!
! (15) ! Ocean					        =>!
! (16) ! Evergreen Shrubs		   	    =>! 5  C4 GRASS
! (17) ! Deciduous Shrubs			    =>! 4  C3 GRASS
! (18) ! Mixed Forest/woodland			=>! 2  CONIFEROUS
! (19) ! Interrupted Forest			    =>! 2  CONIFEROUS
! (20) ! Water and Land Mixtures		=>! 4  C3 GRASS

! LVSTRESS vegetation response type to water stress
! (true:defensive false:offensive) (-)


YDAGS%LVSTRESS(1)=.FALSE.
YDAGS%LVSTRESS(2)=.FALSE.
YDAGS%LVSTRESS(3)=.FALSE.
YDAGS%LVSTRESS(4)=.FALSE.
YDAGS%LVSTRESS(5)=.FALSE. 
YDAGS%LVSTRESS(6)=.FALSE.
YDAGS%LVSTRESS(7)=.FALSE.
YDAGS%LVSTRESS(8)=.FALSE.!not defined 
YDAGS%LVSTRESS(9)=.FALSE.
YDAGS%LVSTRESS(10)=.FALSE.
YDAGS%LVSTRESS(11)=.FALSE.
YDAGS%LVSTRESS(12)=.FALSE.!not defined 
YDAGS%LVSTRESS(13)=.FALSE. 
YDAGS%LVSTRESS(14)=.FALSE.!not defined 
YDAGS%LVSTRESS(15)=.FALSE.!not defined 
YDAGS%LVSTRESS(16)=.FALSE.
YDAGS%LVSTRESS(17)=.FALSE.
YDAGS%LVSTRESS(18)=.FALSE.
YDAGS%LVSTRESS(19)=.FALSE. 
YDAGS%LVSTRESS(20)=.FALSE.
YDAGS%LVSTRESS(0)=  YDAGS%LVSTRESS(7)

! RVR0VT Reference respiration [Kg CO2 m-2 s-1]
!Values optimized by vegetation type with FLUXNET Reco observations

DO JVT=1,NVTYPES ! Loop over PFTs
 DO ICO2TYP=1,2 !C3/C4 plant type
  YDAGS%RVR0VT(JVT,ICO2TYP) = PRVR0VT(JVT,ICO2TYP)
 ENDDO
ENDDO
YDAGS%RVR0VT(0,:) = PRVR0VT(7,:)

! RVCH4QVT reference Q10 for CH4 emission temperature dependency
! optimized at LSCE
! reoptimization at LSCE done on the final stage of CAMS41 (July 2021)
YDAGS%RVCH4QVT(0)=0._JPRB
YDAGS%RVCH4QVT(1)=2.337_JPRB
YDAGS%RVCH4QVT(2)=2.337_JPRB
YDAGS%RVCH4QVT(3)=2.337_JPRB
YDAGS%RVCH4QVT(4)=2.337_JPRB
YDAGS%RVCH4QVT(5)=2.337_JPRB
YDAGS%RVCH4QVT(6)=2.337_JPRB
YDAGS%RVCH4QVT(7)=2.337_JPRB
YDAGS%RVCH4QVT(8)=0._JPRB
YDAGS%RVCH4QVT(9)=2.337_JPRB
YDAGS%RVCH4QVT(10)=2.337_JPRB
YDAGS%RVCH4QVT(11)=2.337_JPRB
YDAGS%RVCH4QVT(12)=0._JPRB
YDAGS%RVCH4QVT(13)=2.337_JPRB
YDAGS%RVCH4QVT(14)=0._JPRB
YDAGS%RVCH4QVT(15)=0._JPRB
YDAGS%RVCH4QVT(16)=2.337_JPRB
YDAGS%RVCH4QVT(17)=2.337_JPRB
YDAGS%RVCH4QVT(18)=2.337_JPRB
YDAGS%RVCH4QVT(19)=2.337_JPRB
YDAGS%RVCH4QVT(20)=0._JPRB

! RVLAIMIN minimum lai (leaf area index) (m2/m2)
! wordt in ISBA uit de ecoclimap LUT gehaald (MINVAL....)
YDAGS%RVLAIMIN(1)=0.3_JPRB
YDAGS%RVLAIMIN(2)=0.3_JPRB
YDAGS%RVLAIMIN(3)=1.0_JPRB
YDAGS%RVLAIMIN(4)=1.0_JPRB
YDAGS%RVLAIMIN(5)=1.0_JPRB
YDAGS%RVLAIMIN(6)=1.0_JPRB
YDAGS%RVLAIMIN(7)=0.3_JPRB
YDAGS%RVLAIMIN(8)=0.0001_JPRB !not defined 
YDAGS%RVLAIMIN(9)=0.3_JPRB
YDAGS%RVLAIMIN(10)=0.3_JPRB
YDAGS%RVLAIMIN(11)=0.3_JPRB
YDAGS%RVLAIMIN(12)=0.0001_JPRB!not defined 
YDAGS%RVLAIMIN(13)=0.3_JPRB
YDAGS%RVLAIMIN(14)=0.0001_JPRB!not defined 
YDAGS%RVLAIMIN(15)=0.0001_JPRB!not defined 
YDAGS%RVLAIMIN(16)=1.0_JPRB
YDAGS%RVLAIMIN(17)=0.3_JPRB
YDAGS%RVLAIMIN(18)=1.0_JPRB
YDAGS%RVLAIMIN(19)=1.0_JPRB
YDAGS%RVLAIMIN(20)=0.3_JPRB
YDAGS%RVLAIMIN(0)=YDAGS%RVLAIMIN(7)

! RVSEFOLD e-folding time for senescence

YDAGS%RVSEFOLD(1)=150._JPRB   ! 150._JPRB*RDAY
YDAGS%RVSEFOLD(2)=150._JPRB   ! 150._JPRB*RDAY
YDAGS%RVSEFOLD(3)=365._JPRB   ! 365._JPRB*RDAY
YDAGS%RVSEFOLD(4)=230._JPRB   !  230._JPRB*RDAY
YDAGS%RVSEFOLD(5)=230._JPRB   ! 230._JPRB*RDAY
YDAGS%RVSEFOLD(6)=365._JPRB   ! 365._JPRB*RDAY
YDAGS%RVSEFOLD(7)=150._JPRB   !  150._JPRB*RDAY
YDAGS%RVSEFOLD(8)=0.001_JPRB  ! 0.001_JPRB*RDAY!not defined 
YDAGS%RVSEFOLD(9)=150._JPRB   ! 150._JPRB*RDAY
YDAGS%RVSEFOLD(10)=150._JPRB  !  150._JPRB*RDAY
YDAGS%RVSEFOLD(11)=150._JPRB  !  150._JPRB*RDAY
YDAGS%RVSEFOLD(12)= 0.001_JPRB!  0.001_JPRB*RDAY!not defined 
YDAGS%RVSEFOLD(13)=150._JPRB  ! 150._JPRB*RDAY
YDAGS%RVSEFOLD(14)=0.001_JPRB ! 0.001_JPRB*RDAY!not defined 
YDAGS%RVSEFOLD(15)=0.001_JPRB ! 0.001_JPRB*RDAY!not defined 
YDAGS%RVSEFOLD(16)=365._JPRB  ! 365._JPRB*RDAY
YDAGS%RVSEFOLD(17)=150._JPRB  ! 150._JPRB*RDAY
YDAGS%RVSEFOLD(18)=365._JPRB  ! 365._JPRB*RDAY
YDAGS%RVSEFOLD(19)=365._JPRB  ! 365._JPRB*RDAY
YDAGS%RVSEFOLD(20)=150._JPRB  ! 150._JPRB*RDAY
YDAGS%RVSEFOLD(0)=YDAGS%RVSEFOLD(7)

! RVGMES mesophyll conductance (m s-1)

! optimized on real biomes with all offensive
YDAGS%RVGMES(1)=.00130_JPRB
YDAGS%RVGMES(2)=.00130_JPRB
YDAGS%RVGMES(3)=.00080_JPRB
YDAGS%RVGMES(4)=.00080_JPRB 
YDAGS%RVGMES(5)=.00140_JPRB
YDAGS%RVGMES(6)=.00110_JPRB  
YDAGS%RVGMES(7)=.00230_JPRB
YDAGS%RVGMES(8)=.000001_JPRB!not defined
YDAGS%RVGMES(9)=.00200_JPRB
YDAGS%RVGMES(10)=.00040_JPRB
YDAGS%RVGMES(11)=.00100_JPRB 
YDAGS%RVGMES(12)=.000001_JPRB!not defined
YDAGS%RVGMES(13)=.00050_JPRB 
YDAGS%RVGMES(14)=.000001_JPRB!not defined
YDAGS%RVGMES(15)=.000001_JPRB!not defined
YDAGS%RVGMES(16)=.00090_JPRB
YDAGS%RVGMES(17)=.00190_JPRB
YDAGS%RVGMES(18)=.00100_JPRB
YDAGS%RVGMES(19)=.00080_JPRB
YDAGS%RVGMES(20)=.00100_JPRB 
YDAGS%RVGMES(0)=YDAGS%RVGMES(7)


! RVGC cuticular conductance (m s-1)
! optimized on real biomes with all offensive
YDAGS%RVGC(1)=0.00015_JPRB
YDAGS%RVGC(2)=0.00020_JPRB
YDAGS%RVGC(3)=0.00020_JPRB
YDAGS%RVGC(4)=0.00020_JPRB 
YDAGS%RVGC(5)=0.0_JPRB
YDAGS%RVGC(6)=0.00025_JPRB
YDAGS%RVGC(7)=0.00020_JPRB
YDAGS%RVGC(8)=0.0000001_JPRB!not defined
YDAGS%RVGC(9)=0.00025_JPRB
YDAGS%RVGC(10)=0.00025_JPRB
YDAGS%RVGC(11)=0.00025_JPRB 
YDAGS%RVGC(12)=0.0000001_JPRB!not defined
YDAGS%RVGC(13)=0.00025_JPRB
YDAGS%RVGC(14)=0.000_JPRB!not defined
YDAGS%RVGC(15)=0.00000001_JPRB!not defined
YDAGS%RVGC(16)=0.00015_JPRB
YDAGS%RVGC(17)=0.00020_JPRB
YDAGS%RVGC(18)=0.0_JPRB 
YDAGS%RVGC(19)=0.00010_JPRB
YDAGS%RVGC(20)=0.00025_JPRB
YDAGS%RVGC(0)=YDAGS%RVGC(7)


! RVDMAX maximum specific humidity deficit tolerated by vegetation (kg kg-1)

YDAGS%RVDMAX(1)=0.050_JPRB
YDAGS%RVDMAX(2)=0.050_JPRB
YDAGS%RVDMAX(3)=0.124_JPRB
YDAGS%RVDMAX(4)=0.124_JPRB 
YDAGS%RVDMAX(5)=0.109_JPRB
YDAGS%RVDMAX(6)=0.124_JPRB
YDAGS%RVDMAX(7)=0.052_JPRB
YDAGS%RVDMAX(8)=0.00001_JPRB!not defined
YDAGS%RVDMAX(9)=0.050_JPRB
YDAGS%RVDMAX(10)=0.050_JPRB
YDAGS%RVDMAX(11)=0.050_JPRB 
YDAGS%RVDMAX(12)=0.00001_JPRB!not defined
YDAGS%RVDMAX(13)=0.050_JPRB
YDAGS%RVDMAX(14)=0.00001_JPRB!not defined
YDAGS%RVDMAX(15)=0.00001_JPRB!not defined
YDAGS%RVDMAX(16)=0.052_JPRB
YDAGS%RVDMAX(17)=0.050_JPRB
YDAGS%RVDMAX(18)=0.124_JPRB 
YDAGS%RVDMAX(19)=0.124_JPRB
YDAGS%RVDMAX(20)=0.050_JPRB
YDAGS%RVDMAX(0)=YDAGS%RVDMAX(7)

! RVF2I critical normilized soil water content for stress parameterisation (-)

DO JVT=0,NVTYPES
  YDAGS%RVF2I(JVT)=0.3_JPRB
ENDDO

! RVCE  specific leaf area (SLA) sensitivity to nitrogen concentration (m2/kg/%)

YDAGS%RVCE(1)=3.79_JPRB
YDAGS%RVCE(2)=5.56_JPRB
YDAGS%RVCE(3)=4.85_JPRB
YDAGS%RVCE(4)=4.85_JPRB 
YDAGS%RVCE(5)=4.83_JPRB
YDAGS%RVCE(6)=4.83_JPRB
YDAGS%RVCE(7)=7.68_JPRB
YDAGS%RVCE(8)=0.0001_JPRB!not defined
YDAGS%RVCE(9)=5.56_JPRB
YDAGS%RVCE(10)=3.79_JPRB
YDAGS%RVCE(11)=5.56_JPRB 
YDAGS%RVCE(12)=0.0001_JPRB!not defined
YDAGS%RVCE(13)=5.56_JPRB
YDAGS%RVCE(14)=0.0001_JPRB!not defined
YDAGS%RVCE(15)=0.0001_JPRB!not defined
YDAGS%RVCE(16)=7.68_JPRB
YDAGS%RVCE(17)=5.56_JPRB
YDAGS%RVCE(18)=4.85_JPRB 
YDAGS%RVCE(19)=4.85_JPRB
YDAGS%RVCE(20)=5.56_JPRB
YDAGS%RVCE(0)=YDAGS%RVCE(7)

! RVCF  lethal minimum value of SLA (m2/kg)

YDAGS%RVCF(1)=9.84_JPRB
YDAGS%RVCF(2)=6.73_JPRB
YDAGS%RVCF(3)=-0.24_JPRB
YDAGS%RVCF(4)=-0.24_JPRB 
YDAGS%RVCF(5)=2.53_JPRB
YDAGS%RVCF(6)=2.53_JPRB
YDAGS%RVCF(7)=-4.33_JPRB
YDAGS%RVCF(8)=0.00001_JPRB!not defined
YDAGS%RVCF(9)=6.73_JPRB
YDAGS%RVCF(10)=9.84_JPRB
YDAGS%RVCF(11)=6.73_JPRB 
YDAGS%RVCF(12)=0.00001_JPRB!not defined
YDAGS%RVCF(13)=6.73_JPRB
YDAGS%RVCF(14)=0.00001_JPRB!not defined
YDAGS%RVCF(15)=0.00001_JPRB!not defined
YDAGS%RVCF(16)=-4.33_JPRB
YDAGS%RVCF(17)=6.73_JPRB
YDAGS%RVCF(18)=-0.24_JPRB 
YDAGS%RVCF(19)=-0.24_JPRB
YDAGS%RVCF(20)=6.73_JPRB
YDAGS%RVCF(0)=YDAGS%RVCF(7)

! RVCNA nitrogen concentration of active biomass (= leaf biomass) (%)

YDAGS%RVCNA(1)=1.3_JPRB
YDAGS%RVCNA(2)=1.3_JPRB 
YDAGS%RVCNA(3)=2.8_JPRB
YDAGS%RVCNA(4)=2.8_JPRB 
YDAGS%RVCNA(5)=2.0_JPRB
YDAGS%RVCNA(6)=2.5_JPRB
YDAGS%RVCNA(7)=1.3_JPRB
YDAGS%RVCNA(8)=0.00001_JPRB!not defined
YDAGS%RVCNA(9)=1.3_JPRB 
YDAGS%RVCNA(10)=1.3_JPRB
YDAGS%RVCNA(11)=1.3_JPRB 
YDAGS%RVCNA(12)=0.00001_JPRB!not defined
YDAGS%RVCNA(13)=1.3_JPRB
YDAGS%RVCNA(14)=0.00001_JPRB!not defined
YDAGS%RVCNA(15)=0.00001_JPRB!not defined
YDAGS%RVCNA(16)=1.3_JPRB 
YDAGS%RVCNA(17)=1.3_JPRB
YDAGS%RVCNA(18)=2.8_JPRB 
YDAGS%RVCNA(19)=2.8_JPRB
YDAGS%RVCNA(20)=1.3_JPRB
YDAGS%RVCNA(0)=YDAGS%RVCNA(7)

! RVBSLAI ratio d(biomass)/d(lai) (kg/m2) 

YDAGS%RVBSLAI(1)=.06_JPRB
YDAGS%RVBSLAI(2)=.08_JPRB
YDAGS%RVBSLAI(3)=.50_JPRB
YDAGS%RVBSLAI(4)=.50_JPRB
YDAGS%RVBSLAI(5)=.125_JPRB
YDAGS%RVBSLAI(6)=.25_JPRB
YDAGS%RVBSLAI(7)=.08_JPRB
YDAGS%RVBSLAI(8)=.00001_JPRB!not defined
YDAGS%RVBSLAI(9)=.08_JPRB
YDAGS%RVBSLAI(10)=.06_JPRB
YDAGS%RVBSLAI(11)=.08_JPRB
YDAGS%RVBSLAI(12)=.00001_JPRB!not defined
YDAGS%RVBSLAI(13)=.08_JPRB
YDAGS%RVBSLAI(14)=.00001_JPRB!not defined
YDAGS%RVBSLAI(15)=.00001_JPRB!not defined
YDAGS%RVBSLAI(16)=.08_JPRB
YDAGS%RVBSLAI(17)=.08_JPRB
YDAGS%RVBSLAI(18)=.5_JPRB
YDAGS%RVBSLAI(19)=.5_JPRB
YDAGS%RVBSLAI(20)=.08_JPRB
YDAGS%RVBSLAI(0)=YDAGS%RVBSLAI(7)

DO JVT=0,NVTYPES
  IF (JVT==3 .OR. JVT==4 .OR. JVT==5 .OR. JVT==6 .OR. JVT==18 .OR. JVT==19) THEN  !BATS classification matching
    YDAGS%RXBOMEGAM(JVT)=RXBOMEGA*1.0_JPRB
  ELSE
    YDAGS%RXBOMEGAM(JVT)=RXBOMEGA

  ENDIF
ENDDO
!    -------------------------------------------------------------------------

!          1.3 INITIALIZE PARAMETERS DEPENDING ON PHOTOSYNTHESIS MECHANISM
!              ---------- ---------- --------- -- -------------- ---------
IF (.NOT.ALLOCATED(YDAGS%RVTOPT)   ) ALLOCATE (YDAGS%RVTOPT   (0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVFZERO)  ) ALLOCATE (YDAGS%RVFZERO  (0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVFZEROST)) ALLOCATE (YDAGS%RVFZEROST(0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVEPSO)   ) ALLOCATE (YDAGS%RVEPSO   (0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVGAMM)   ) ALLOCATE (YDAGS%RVGAMM   (0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVQDGAMM) ) ALLOCATE (YDAGS%RVQDGAMM (0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVQDGMES) ) ALLOCATE (YDAGS%RVQDGMES (0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVT1GMES) ) ALLOCATE (YDAGS%RVT1GMES (0:NVTYPES))
IF (.NOT.ALLOCATED(YDAGS%RVT2GMES) ) ALLOCATE (YDAGS%RVT2GMES (0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVAMMAX)  ) ALLOCATE (YDAGS%RVAMMAX  (0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVQDAMMAX)) ALLOCATE (YDAGS%RVQDAMMAX(0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVT1AMMAX)) ALLOCATE (YDAGS%RVT1AMMAX(0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVT2AMMAX)) ALLOCATE (YDAGS%RVT2AMMAX(0:NVTYPES)) 
IF (.NOT.ALLOCATED(YDAGS%RVAH)     ) ALLOCATE (YDAGS%RVAH     (0:NVTYPES))
IF (.NOT.ALLOCATED(YDAGS%RVBH)     ) ALLOCATE (YDAGS%RVBH     (0:NVTYPES)) 

DO JVT=0,NVTYPES
  IF (JVT==7 .OR. JVT==16 .OR. JVT==0) THEN  !BATS classification matching
    ICO2TYP = 2   ! C4 type
  ELSE
    ICO2TYP = 1   ! C3 type
  END IF
  YDAGS%RVTOPT    (JVT) = ZTOPT    (ICO2TYP)
  YDAGS%RVFZERO   (JVT) = ZFZERO   (ICO2TYP)
  YDAGS%RVFZEROST (JVT) = ZFZEROST (ICO2TYP)  
  YDAGS%RVEPSO    (JVT) = ZEPSO    (ICO2TYP)
  YDAGS%RVGAMM    (JVT) = ZGAMM    (ICO2TYP)
  YDAGS%RVQDGAMM  (JVT) = ZQDGAMM  (ICO2TYP)
  YDAGS%RVQDGMES  (JVT) = ZQDGMES  (ICO2TYP)
  YDAGS%RVT1GMES  (JVT) = ZT1GMES  (ICO2TYP)
  YDAGS%RVT2GMES  (JVT) = ZT2GMES  (ICO2TYP)
  YDAGS%RVAMMAX   (JVT) = ZAMMAX   (ICO2TYP)
  YDAGS%RVQDAMMAX (JVT) = ZQDAMMAX (ICO2TYP)
  YDAGS%RVT1AMMAX (JVT) = ZT1AMMAX (ICO2TYP)
  YDAGS%RVT2AMMAX (JVT) = ZT2AMMAX (ICO2TYP)
  YDAGS%RVAH      (JVT) = ZAH      (ICO2TYP)
  YDAGS%RVBH      (JVT) = ZBH      (ICO2TYP)
ENDDO

! optimized on real biomes with all offensive
 YDAGS%RVAMMAX(1)=2.2E-6_JPRB
 YDAGS%RVAMMAX(2)=3.0E-6_JPRB
 YDAGS%RVAMMAX(3)=2.2E-6_JPRB
 YDAGS%RVAMMAX(4)=2.2E-6_JPRB
 YDAGS%RVAMMAX(5)=1.833E-6_JPRB
 YDAGS%RVAMMAX(6)=1.833E-6_JPRB
 YDAGS%RVAMMAX(7)=1.833E-6_JPRB
 YDAGS%RVAMMAX(8)=1.833E-6_JPRB
 YDAGS%RVAMMAX(9)=3.0E-6_JPRB
 YDAGS%RVAMMAX(10)=1.833E-6_JPRB
 YDAGS%RVAMMAX(11)=1.833E-6_JPRB
 YDAGS%RVAMMAX(12)=1.833E-6_JPRB
 YDAGS%RVAMMAX(13)=1.833E-6_JPRB
 YDAGS%RVAMMAX(14)=1.833E-6_JPRB
 YDAGS%RVAMMAX(15)=1.833E-6_JPRB
 YDAGS%RVAMMAX(16)=1.833E-6_JPRB
 YDAGS%RVAMMAX(17)=1.833E-6_JPRB
 YDAGS%RVAMMAX(18)=2.2E-6_JPRB
 YDAGS%RVAMMAX(19)=2.2E-6_JPRB
 YDAGS%RVAMMAX(20)=1.833E-6_JPRB
 YDAGS%RVAMMAX(0)= YDAGS%RVAMMAX(7)


! optimized on real biomes with all offensive
 YDAGS%RVFZEROST(1)=0.85_JPRB
 YDAGS%RVFZEROST(2)=0.65_JPRB
 YDAGS%RVFZEROST(3)=0.95_JPRB
 YDAGS%RVFZEROST(4)=0.95_JPRB
 YDAGS%RVFZEROST(5)=0.95_JPRB
 YDAGS%RVFZEROST(6)=0.95_JPRB
 YDAGS%RVFZEROST(7)=0.70_JPRB
 YDAGS%RVFZEROST(8)=0.95_JPRB
 YDAGS%RVFZEROST(9)=0.95_JPRB
 YDAGS%RVFZEROST(10)=0.92_JPRB
 YDAGS%RVFZEROST(11)=0.80_JPRB
 YDAGS%RVFZEROST(12)=0.95_JPRB
 YDAGS%RVFZEROST(13)=0.96_JPRB
 YDAGS%RVFZEROST(14)=0.95_JPRB
 YDAGS%RVFZEROST(15)=0.95_JPRB
 YDAGS%RVFZEROST(16)=0.72_JPRB
 YDAGS%RVFZEROST(17)=0.96_JPRB
 YDAGS%RVFZEROST(18)=0.95_JPRB
 YDAGS%RVFZEROST(19)=0.95_JPRB
 YDAGS%RVFZEROST(20)=0.95_JPRB
 YDAGS%RVFZEROST(0)= YDAGS%RVFZEROST(7)


!    Calculates the Gaussian weights for integration of net assimilation
!    and stomatal conductance over the canopy depth
!    using Gaussian Quadrature
!     *KN*        NUMBER OF POINTS AT WHICH GAUSSIAN WEIGHTS       -
!                 ARE EVALUATED/NEEDED    
!     *ZX1*       MATHEMATICAL/NUMERICAL VALUES NEEDED FOR         -
!     *ZX2*                             WEIGHT COMPUTATION
!     OUTPUT PARAMETERS (REAL):
!     *ZX*        ABSCISSA                                         -
!     *ZW*        GAUSSIAN WIEGHTS                                 -
!--------------
ZX1=0.0_JPRB
ZX2=1.0_JPRB
KN=SIZE(RABC)
!--------------  
IM=(KN+1)/2
ZXM=0.5_JPRB*(ZX2+ZX1)
ZXL=0.5_JPRB*(ZX2-ZX1)
         
DO JI=1,IM
      ZZ=COS(RPI*(FLOAT(JI)-.25_JPRB)/(FLOAT(KN)+.5_JPRB))
   !  begin iteration:
      DO JK=1,100
         ZP1=1._JPRB
         ZP2=0._JPRB
         DO JJ=1,KN
            ZP3=ZP2
            ZP2=ZP1
            ZP1=((2._JPRB*(JJ)-1._JPRB)*ZZ*ZP2-(FLOAT(JJ)-1._JPRB)*ZP3)/JJ
         END DO
         ZPP=FLOAT(KN)*(ZZ*ZP1-ZP2)/(ZZ*ZZ-1._JPRB)
         ZZ1=ZZ
         ZZ=ZZ1-ZP1/ZPP
         IF(ABS(ZZ-ZZ1).LE.ZPEPS)EXIT
      ENDDO 
   !  end iteration.
   !  compute abscissa
       RABC(JI)=ZXM-ZXL*ZZ
       RABC(KN+1-JI)=ZXM+ZXL*ZZ   
   !  compute weights
      RPOI(JI)=2._JPRB*ZXL/((1._JPRB-ZZ*ZZ)*ZPP*ZPP)
      RPOI(KN+1-JI)=RPOI(JI)

ENDDO 
!-------------- 

END ASSOCIATE
IF (LHOOK) CALL DR_HOOK('SUCOTWO_MOD:SUCOTWO',1,ZHOOK_HANDLE)

END SUBROUTINE SUCOTWO
END MODULE SUCOTWO_MOD
 
